Name, Vorname:
Hofmann Fabian
Richter Daniel
2. Schätzung der Fourierreihe durch FFT
3. Schätzung des Fourierspektrums von Sinustönen mit überlagertem Rauschen durch FFT
4. Fourierspektrum eines Rechteckpulses
Nach dem Praktikumsversuch exportieren Sie das Notebook mit Textantworten, Codezellen und Plots als HTML (File -> Export Notebook As ... -> Export Notebook to HTML) und reichen es in Moodle ein.</div>
import os, sys
module_path = os.path.abspath(os.path.join('..')) # append directory one level up to import path
if module_path not in sys.path: # ... if it hasn't been appended already
sys.path.append(module_path)
import dsp_fpga_lib as dsp
dsp.versions() # print versions
%matplotlib inline
import matplotlib.pyplot as plt
size = {"figsize":(12,5)} # Plotgröße in Inch
import numpy as np
import scipy.signal as sig
import wave
from IPython.display import Audio, display
#-----------------------------------------------------------------------------
def wav2np(filename):
""" Read the wav-file and convert it to a one or two-dimensional numpy array,
depending on the number of channels.
Properties of the WAV-file are stored as function attributes (evil)
"""
wf = wave.open(filename,'rb')
wav2np.N_CH = wf.getnchannels() # number of channels
wav2np.W = wf.getsampwidth() # wordlength per sample in bytes
wav2np.N_FR = wf.getnframes() # number of frames
wav2np.f_S = wf.getframerate() # sample (frame) rate
print("{0} channels with {1} frames of {2} bytes and f_S = {3} Hz.".format(wav2np.N_CH, wav2np.N_FR, wav2np.W, wav2np.f_S))
if wav2np.W == 2:
samples_in = np.frombuffer(wf.readframes(-1), dtype=np.int16) # read wav data as 16 bit integers, R and L samples are interleaved
elif wav2np.W == 1:
samples_in = np.frombuffer(wf.readframes(-1), dtype=np.int8) # read wav data as 8 bit integers, R and L samples are interleaved
else:
raise TypeError("Unknown data format: {0} bytes".format(wav2np.W))
samples = np.array([samples_in[idx::wav2np.N_CH] for idx in range(wav2np.N_CH)], dtype=np.int32) # deinterleave channels to numpy array N_CHAN x N_FRAMES
return samples
Python version: 3.10.8 Numpy: 1.23.5 Scipy: 1.9.3 Matplotlib: 3.6.2 module://matplotlib_inline.backend_inline
In diesem Abschnitt analysieren Sie die beiden Filter aus der folgenden Abbildung:
| Filter 1 | Filter 2 |
Für die Simulation müssen wir zuerst die Filter als Koeffizientenvektoren a und b darstellen (siehe folgendes Bild und Abschnitt 3 von LAB 1).
Als "Codesteinbruch" können Sie wieder das Notebook 02_LTF/LTF-Filter_properties.ipynb verwenden.
Alternativ verwenden Sie pyfda: Im Tab "b,a" importieren Sie Koeffizienten aus einem CSV-File oder (nach Auswahl von "Clipboard" in den Einstellungen) direkt aus der Zwischenablage. Die Koeffizienten b des nicht-rekursiven Teils geben Sie hierfür getrennt durch Kommata an, optional in einer zweiten Zeile die Koeffizienten a des rekursiven Teils.
Stellen Sie Filter 1 und Filter 2 (s.o.) als Koeffizientenvektoren a und b dar. Informationen dazu finden Sie im Abschnitt 3 von LAB 1 und in obiger Abbildung. Beachten Sie:
1j, dementsprechend werden imaginäre Zahlen z.B. als 0.3j repräsentiert.Testen Sie mit dem P/N Diagramm, ob Rechnung und Simulation zusammen passen.
Lassen Sie sich die Impulsantwort anzeigen, für das FIR-Filter ist das relativ einfach ;-), beim IIR-Filter benötigen Sie dsp.impz()
Plotten Sie Betragsfrequenzgang, Phasengang und Gruppenlaufzeit.
#Filter
a_FIR = [1, 0, 0, 0]
b_FIR = [0, 0.5, 1, 0.5]
#Filter
a_IIR = [1, 0, 0.9]
b_IIR = [0, 0, 1]
#---------------
fig, (ax1,ax2) = plt.subplots(nrows=1, ncols=2,**size)
ax1.set_title('P/N Diagramm: FIR')
dsp.zplane(b_FIR,a_FIR, plt_ax=ax1);
print("Nullstellen FIR: {0}".format(np.roots(b_FIR)))
ax2.set_title('P/N Diagramm: IIR')
dsp.zplane(b_IIR,a_IIR, plt_ax=ax2);
print("Nullstellen IIR: {0}".format(np.roots(b_IIR)))
if type(a_IIR) in {list, np.ndarray} and len(a_IIR) > 1:
print("Polstellen IIR: {0}\n".format(np.roots(a_IIR)))
Nullstellen FIR: [-1. -1.] Nullstellen IIR: [] Polstellen IIR: [-0.+0.9486833j 0.-0.9486833j]
Der Befehl h,t = dsp.impz(b,a) berechnet die Impulsantwort, die Sie dann plotten können, am Besten als stem-Plot.
Warum ist beim IIR-Filter jeder zweite Impuls Null?
fig, (ax1,ax2) = plt.subplots(nrows=2, ncols=1,**size)
ax1.set_title('Impulsantwort FIR')
h_FIR,t_FIR = dsp.impz(b_FIR,a_FIR);
ax1.stem(t_FIR, h_FIR)
ax2.set_title('Impulsantwort IIR')
h_IIR,t_IIR = dsp.impz(b_IIR,a_IIR);
ax2.stem(t_IIR, h_IIR)
fig.set_tight_layout(True)
Aus dem komplexen Frequenzgang omega, H = sig.freqz(b,a) ermitteln Sie mit np.abs() und np.angle() Betrags- und Phasengang. Defaultmäßig werden 512 Frequenzpunkte zwischen 0 und $f_S/2$ bestimmt und als normierte Kreisfrequenz zurückgegeben.
Die Gruppenlaufzeit ermitteln Sie mit w, H = sig.group_delay((b,a), omega). Mit omega übergeben Sie die Frequenzen, an denen die Gruppenlaufzeit berechnet werden soll (z.B. die, die Sie aus der Berechnung des komplexen Frequenzgangs erhalten haben).
Der Plot der Gruppenlaufzeit sieht u.U. etwas seltsam aus, mit set_ylim([y_min, y_max]) passen Sie die Grenzen an.
fig, ax = plt.subplots(nrows=3, ncols=2, figsize=(12,8))
#FIR
w_FIR, H_FIR = sig.freqz(b_FIR, a_FIR)
wgr_FIR, tau_FIR = sig.group_delay((b_FIR, a_FIR), w_FIR)
#IIR
w_IIR, H_IIR = sig.freqz(b_IIR, a_IIR)
wgr_IIR, tau_IIR = sig.group_delay((b_IIR, a_IIR), w_IIR)
#FIR
#Betrag
ax[0][0].set_title("FIR-Betrag")
ax[0][0].plot(w_FIR / (2*np.pi), np.abs(H_FIR))
ax[0][0].set_ylabel(r"$|H| \;\rightarrow$")
ax[0][0].set_xlabel(r"$F \;\rightarrow$")
#Phase
ax[1][0].set_title("FIR-Phase")
ax[1][0].plot(w_FIR / (2*np.pi), np.angle(H_FIR))
ax[1][0].set_ylabel(r"$Ph° \;\rightarrow$")
ax[1][0].set_xlabel(r"$F \;\rightarrow$")
#Gruppenlaufzeit
ax[2][0].set_title("FIR-Gruppenlaufzeit")
ax[2][0].plot(wgr_FIR / (2*np.pi), tau_FIR)
ax[2][0].set_ylabel(r"$T_g \;\rightarrow$")
ax[2][0].set_xlabel(r"$F \;\rightarrow$")
#IIR
ax[0][1].set_title("IIR")
ax[0][1].plot(w_IIR / (2*np.pi), np.abs(H_IIR))
ax[0][1].set_ylabel(r"$|H| \;\rightarrow$")
ax[0][1].set_xlabel(r"$F \;\rightarrow$")
#Phase
ax[1][1].set_title("IIR-Phase")
ax[1][1].plot(w_IIR / (2*np.pi), np.angle(H_IIR))
ax[1][1].set_ylabel(r"$Ph° \;\rightarrow$")
ax[1][1].set_xlabel(r"$F \;\rightarrow$")
#Gruppenlaufzeit
ax[2][1].set_title("IIR-Gruppenlaufzeit")
ax[2][1].plot(wgr_IIR / (2*np.pi), tau_IIR)
ax[2][1].set_ylabel(r"$T_g \;\rightarrow$")
ax[2][1].set_xlabel(r"$F \;\rightarrow$")
fig.set_tight_layout(True)
#set_ylim([y_min, y_max])
Im folgenden filtern wir Sprach- oder Rauschsignale mit unseren beiden Filtern und hören und schauen uns das Resultat an. Eine FIR-Filterung könnten Sie wieder mit convolve(x,h) durchführen (nur für eindimensionale Arrays). Warum funktioniert bei IIR-Filtern die convolve-Methode nicht?
Für IIR und FIR-Filter können Sie die Routine y = sig.lfilter(b,a,x) verwenden (funktioniert auch mit zweidimensionalen Arrays).
Filtern Sie ein Rauschsignal oder einen WAV-File aus dem Unterordner medien mit dem FIR- und dem IIR-Filter. Rauschen erzeugen Sie wieder z.B. mit x_n = np.random.randn(16000), WAV-Files wandeln Sie mit der Hilfsfunktion wav2array(filename) in ein ein- oder zweidimensionales Array um.
Hören Sie sich die Wirkung der beiden Filter an mit der Audio Klasse aus dem IPython.display - Modul:
display(Audio((data=None, rate=None), data kann dabei ein ein- oder zweidimensionales numpy-Array oder Liste sein, ein Filename oder auch eine URL. Der Parameter rate definiert die Abtastrate.
Schauen Sie sich einen Ausschnitt des ursprünglichen und des gefilterten Signals im gleichen Plot an und vergleichen Sie.
x_n = np.random.randn(16000)
x_w = wav2np("../medien/87778__marcgascon7__vocals.wav")
print(np.shape(x_w))
y_FIR = sig.lfilter(b_FIR,a_FIR,x_w)
y_IIR = sig.lfilter(b_IIR,a_IIR,x_w)
print("Orginal")
display(Audio(data=x_w, rate=wav2np.f_S))
print("FIR Filter")
display(Audio(data=y_FIR, rate=wav2np.f_S))
print("IIR Filter")
display(Audio(data=y_IIR, rate=wav2np.f_S))
2 channels with 349952 frames of 2 bytes and f_S = 44100 Hz. (2, 349952) Orginal
FIR Filter
IIR Filter
fig, ax = plt.subplots(**size)
l = 4000
start = 20000
t = np.arange(start, start + l, 1)
ax.plot(t, x_w[0][start:start + l:1], label="x_w")
ax.plot(t, y_FIR[0][start:start + l:1], label="y_FIR")
ax.plot(t, y_IIR[0][start:start + l:1], label="y_IIR")
ax.legend();
In diesem Versuchsteil bestimmen wir die Fourierkoeffizienten einer rechteckförmigen zeitkontinuierlichen Pulsfolge mit $T_1 = 1\,$ms. Der Duty Cycle $\alpha = 0.25$ und Amplitude $A=1000\,$mV, die Pulsbreite ist also $\Delta T = T_1/4$ (siehe nächster Plot).
fig, ax = plt.subplots(**size)
t = np.arange(-0.2e-3, 4e-3, 1/640e3) # pseudo-analoger Zeitvektor in ms (mit 640 kHz so fein abgetastet, dass der Unterschied nicht auffällt)
y = 500*sig.square(t*1e3*2*np.pi+np.pi/4, duty=0.25) + 500
ax.plot(t, y, 'r', lw=2)
ax.axvline(0, ls='-',color='k')
ax.set_xlabel(r"$t \; / \;\mathrm{ms} \;\rightarrow$")
ax.set_ylabel(r"$y(t)\; / \;\mathrm{mV} \;\rightarrow$");
Ein mit $T_1$ periodisches Signal $y(t)$ lässt sich mit Hilfe der Fourierreihe für reellwertige Signale (Index $k \in \mathbb{N}_0$) darstellen:
$$ y(t)= {a_0} + 2\sum_{k=1}^\infty a_k \cos 2 \pi k f_1 t + b_k \sin 2 \pi k f_1 t \; k \in \mathbb{N} $$Zunächst berechnen wir die Koeffizienten $a_k$, $b_k$ für $k = 0, 1, \ldots, 11$ mit Hilfe der Formel für die Fourierreihe von reellwertigen Zeitsignalen. Da das Signal achsensymmetrisch ist, sind die imaginärwertigen Koeffizienten $b_k=0$.
$$ \begin{align}{a}_k &= \frac 1 {T_1} \int_{-T_1/2}^{T_1/2}y(t)\cos(2 k \pi f_1 t) \, \text{d}t = \frac {A} {T_1} \int_{-\alpha T_1/2}^{\alpha T_1/2}\cos(2 k \pi f_1 t) \, \text{d}t \\ &= \left. \frac {A} {T_1 \cdot 2 k \pi f_1} \sin(2 k \pi f_1 t)\, \right|_{-\alpha T_1/2}^{\alpha T_1/2} = {\frac {A} { k \pi} \sin ( k \alpha \pi) }= { {A \alpha}\, \mathrm{si} ( k \alpha \pi) } \end{align}$$Schreiben Sie ein Skript, das die ersten 11 Koeffizienten der Rechteckfunktion exakt berechnet. Welchen Vorteil hat es, die sinc(x) Funktion zu verwenden anstelle von sin(x)/x? Aber Achtung: In Numpy (und in Matlab) ist die sinc-Funktion definiert als $\mathrm{sinc}(x) = sin(\pi x)/\pi x$, das $\pi$ im Argument müssen Sie daher weglassen. Drucken Sie die Ergebnisse übersichtlich in eine Tabelle aus (s.u.)
Einen formatierten Ausdruck erhalten Sie z.B. mit print("\n{0:7.2f} | ".format(Y[i]), end="") (der Teil in der geschweiften Klammer wird ersetzt durch Y[i] und formatiert mit insgesamt 7 Stellen, 2 Nachkommastellen, keinen Zeilenumbruch). Eingebettet in eine for Schleife bekommen Sie so schnell eine übersichtliche Tabelle.
print(a,end="") ersetzt den Zeilenumbruch am Ende des Printbefehls durch ein anderes Zeichen (oder nichts wie hier)"{0} ich heiße und bin {1} Jahre alt".format("Yoda", 877).Der Ausdruck in der geschweiften Klammer kann durch Formatierungsanweisung ergänzt werden, {0:7.2f} formatiert das Argument 0 als Float mit 2 Nachkommastellen und insgesamt mindestens 7 Stellen. {0:>7} lässt ebenfalls Platz für mindestens 7 Stellen oder Zeichen und richtet den Ausdruck rechtsbündig aus.
Auf https://www.python-kurs.eu/python3_formatierte_ausgabe.php finden Sie eine übersichtliche Grafik hierzu.
k=np.arange(16)
alpha = 4/16; A = 1000
Y=[]
print("i = | ", end="")
for i in k:
Y.append(A*alpha*np.sinc(k[i]*alpha)) #*np.pi)
print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for i in k:
print("{0:>7.2f} | ".format(Y[i]), end="")
i = | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | Y = | 250.00 | 225.08 | 159.15 | 75.03 | 0.00 | -45.02 | -53.05 | -32.15 | -0.00 | 25.01 | 31.83 | 20.46 | 0.00 | -17.31 | -22.74 | -15.01 |
Um das Spektrum des zeitkontinuierlichen Signals mit einer FFT abschätzen zu können, muss es zunächst abgetastet werden und zwar über eine ganzzahlige Anzahl Perioden $T_1$ (da es ja periodisch ist), wir starten zunächst mit einer Periode, $T_{mess1} = T_1$. Für eine effiziente Berechnung sollen $N_1 = 2^4 = 16$ Samples genommen werden.
Welche Abtastfrequenz $f_{S1}$ ist dafür erforderlich?
Zunächst einmal benötigen Sie einen passenden Zeitvektor z.B. mit t1 = np.arange(N_1)/f_S1 oder durch Abtasten des "zeitkontinuierlichen" Zeitvektors t (s.u.). Das abgetastete Signal $y_1[n]$ erzeugen Sie:
Händisch mit y1 = np.ones(N_1); y1[4:N_1] = 0 oder durch y1 = np.concatenate((np.ones(4), np.zeros(12)))
Durch Abtasten des "zeitkontinuierlichen" Signals ($f_S = 640 \text{ kHz}$) mit y1 = y[::40], mit diesem Befehl wird von y nur jedes vierzigste Element nach y1 kopiert. Das funktioniert natürlich nur, wenn wie hier $f_{S1} = f_S / 40$ ist. Wenn wie hier Start- und Endwert fehlen, werden alle Elemente von y vom ersten bis zum letzten berücksichtigt, Sie müssen ggf. noch Start- und Endwert anpassen, um den richtigen Ausschnitt finden.
Durch Berechnung aus dem Zeitvektor mit Hilfe der np.where() Funktion wie im vorigen Lab.
Egal für welche Variante Sie sich entschieden haben, plotten Sie das abgetastete Signal mit dem folgenden Skript.
fig, ax = plt.subplots(**size)
N_1 = 16
f_S1 = 16e3
fs = 640000
t1 = np.arange(N_1)/f_S1
y1 = y[:N_1*40:40]
#print(np.shape(t1))
#print(np.shape(y1))
ax.stem(t1, y1, 'r',use_line_collection=True)
ax.set_xlabel(r"$k T_S \; / \;\mathrm{s} \;\rightarrow$")
ax.set_ylabel(r"$y(k T_S)\; / \;\mathrm{mV} \;\rightarrow$");
C:\Users\Fabia\AppData\Local\Temp\ipykernel_2544\2715965866.py:9: MatplotlibDeprecationWarning: The 'use_line_collection' parameter of stem() was deprecated in Matplotlib 3.6 and will be removed two minor releases later. If any parameter follows 'use_line_collection', they should be passed as keyword, not positionally. ax.stem(t1, y1, 'r',use_line_collection=True)
In numpy berechnen wir die FFT mit Y = np.fft.fft(y, N) (die FFT liefert komplexe Werte ...). Ohne den Parameter N wird die FFT über das gesamte Signal berechnet mit der Anzahl der Datenpunkte. Nutzen Sie das Notebook 03_DFT/DFT-Skalierung.ipynb zur korrekten Berechnung und Skalierung des Signals. Mit f = np.fft.fftfreq(N, T_S) erzeugen Sie eine passenden Frequenzvektor.
fig, ax = plt.subplots(**size)
N1 = 16
f_S1 = N1 * 1e3
#t1 = np.arange(N1)/f_S1
#t1 = 0
#print(y1)
#FFT
Y1 = np.fft.fft(y1, N1)
Y1 = ((Y1)/N1)[:N1:] #(np.absolute(Y1)/N1)[:N1:]
f1 = np.fft.fftfreq(N1, 1/f_S1)
#Plot
ax.stem(f1, Y1, 'r',use_line_collection=True)
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{mV} \;\rightarrow$");
#ax.stem(f1[:12], Y, 'r',use_line_collection=True)
#ax.set_xlabel(r"$k T_S \; / \;\mathrm{s} \;\rightarrow$")
#ax.set_ylabel(r"$y(k T_S)\; / \;\mathrm{mV} \;\rightarrow$");
# ----------- Tabellenausgabe --------------
Y1 = np.absolute(Y1)
print("FFT:")
print("i = | ", end="")
for i in range(N1):
print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for i in range(N1):
print("{0:>7.2f} | ".format(Y1[i]), end="")
print("\n---------------------------------------------------------------------------------------------------------------------------------------------------------------------\nPer Hand:")
print("i = | ", end="")
for i in k:
print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for i in k:
print("{0:>7.2f} | ".format(Y[i]), end="")
C:\Users\Fabia\AppData\Local\Temp\ipykernel_2544\4006529774.py:14: MatplotlibDeprecationWarning: The 'use_line_collection' parameter of stem() was deprecated in Matplotlib 3.6 and will be removed two minor releases later. If any parameter follows 'use_line_collection', they should be passed as keyword, not positionally. ax.stem(f1, Y1, 'r',use_line_collection=True)
FFT: i = | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | Y = | 250.00 | 226.53 | 163.32 | 79.55 | 0.00 | 53.15 | 67.65 | 45.06 | 0.00 | 45.06 | 67.65 | 53.15 | 0.00 | 79.55 | 163.32 | 226.53 | --------------------------------------------------------------------------------------------------------------------------------------------------------------------- Per Hand: i = | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | Y = | 250.00 | 225.08 | 159.15 | 75.03 | 0.00 | 45.02 | 53.05 | 32.15 | 0.00 | 25.01 | 31.83 | 20.46 | 0.00 | 17.31 | 22.74 | 15.01 |
Die Abtastfrequenz wird jetzt vervierfacht auf den Wert $f_{S2} = 64 \text{ kHz}$.
Kopieren Sie das Skript aus dem vorigen Unterpunkt und passsen Sie es an die ermittelten Parameter an.
Vergleichen Sie erneut Rechnung mit Simulation
fig, ax = plt.subplots(**size)
N2 = 64
f_S2 = N2 * 1e3
y2 = y[:N2*10:10]
Y2 = np.fft.fft(y2, N2)
Y2 = (np.absolute(Y2)/N2)[:N2:]
f2 = np.fft.fftfreq(N2, 1/f_S2)
#Plot
ax.stem(f2, Y2, 'r',use_line_collection=True)
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{mV} \;\rightarrow$");
#ax.stem(f1[:12], Y, 'r',use_line_collection=True)
#ax.set_xlabel(r"$k T_S \; / \;\mathrm{s} \;\rightarrow$")
#ax.set_ylabel(r"$y(k T_S)\; / \;\mathrm{mV} \;\rightarrow$");
# ----------- Tabellenausgabe --------------
print("FFT1:")
print("i = | ", end="")
for i in range(N1):
print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for i in range(N1):
print("{0:>7.2f} | ".format(Y1[i]), end="")
print("\n---------------------------------------------------------------------------------------------------------------------------------------------------------------------\nFFT2:")
#Y2 = np.absolute(Y2)/64
print("i = | ", end="")
for i in range(N1):
print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for i in range(N1):
print("{0:>7.2f} | ".format(Y2[i]), end="")
print("\n---------------------------------------------------------------------------------------------------------------------------------------------------------------------\nPer Hand:")
print("i = | ", end="")
Y = np.absolute(Y)
for i in k:
print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for i in k:
print("{0:>7.2f} | ".format(Y[i]), end="")
FFT1: i = | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | Y = | 250.00 | 226.53 | 163.32 | 79.55 | 0.00 | 53.15 | 67.65 | 45.06 | 0.00 | 45.06 | 67.65 | 53.15 | 0.00 | 79.55 | 163.32 | 226.53 | --------------------------------------------------------------------------------------------------------------------------------------------------------------------- FFT2: i = | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | Y = | 250.00 | 225.17 | 159.41 | 75.30 | 0.00 | 45.47 | 53.83 | 32.80 | 0.00 | 25.84 | 33.15 | 21.49 | 0.00 | 18.55 | 24.63 | 16.45 | --------------------------------------------------------------------------------------------------------------------------------------------------------------------- Per Hand: i = | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | Y = | 250.00 | 225.08 | 159.15 | 75.03 | 0.00 | 45.02 | 53.05 | 32.15 | 0.00 | 25.01 | 31.83 | 20.46 | 0.00 | 17.31 | 22.74 | 15.01 |
C:\Users\Fabia\AppData\Local\Temp\ipykernel_2544\2007472309.py:11: MatplotlibDeprecationWarning: The 'use_line_collection' parameter of stem() was deprecated in Matplotlib 3.6 and will be removed two minor releases later. If any parameter follows 'use_line_collection', they should be passed as keyword, not positionally. ax.stem(f2, Y2, 'r',use_line_collection=True)
Die Abtastfrequenz ist wieder wie zu Beginn $f_{S3} = 16 \text{ kHz}$, jetzt soll getestet werden wie sich eine Verlängerung des Messfensters um den Faktor 4 auswirkt.
Da Ihr Messfenster jetzt mehrere Perioden des Signals umfasst, müssen Sie Ihr Signal erzeugen entweder
Händisch durch mehrfaches Aneinanderhängen eines Pulses z.B. mit y3 = np.tile(y1, 4)
Durch Abtasten des "zeitkontinuierlichen" Signals wie am Anfang von 2.2 beschrieben. Achten Sie wieder darauf, Anfang und Ende passend zu wählen.
Plotten Sie Ihr Zeitsignal, damit Sie sicher sind dass Ihr Signal so aussieht wie Sie sich das vorstellen.
fig, ax = plt.subplots(**size)
N3 = 64
f_S3 = 16 * 1e3
t3 = np.arange(N3)/f_S3
y3 = y[:N3*40:40]
#Plot Signal gefenstert
ax.stem(t3, y3, 'r',use_line_collection=True)
ax.set_xlabel(r"$t\; / \;\mathrm{s} \;\rightarrow$")
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{mV} \;\rightarrow$");
#FFT
Y3 = np.fft.fft(y3, N3)
Y3 = (np.absolute(Y3)/N3)
f3 = np.fft.fftfreq(N3, 1/f_S3)
#Plot FFT
fig, ax = plt.subplots(**size)
ax.stem(f3, Y3, 'r',use_line_collection=True)
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{mV} \;\rightarrow$");
#fig, ax = plt.subplots(**size)
#ax.stem(f1[:12], Y, 'r',use_line_collection=True)
#ax.set_xlabel(r"$k T_S \; / \;\mathrm{s} \;\rightarrow$")
#ax.set_ylabel(r"$y(k T_S)\; / \;\mathrm{mV} \;\rightarrow$");
# ----------- Tabellenausgabe --------------
print("FFT1:")
print("i = | ", end="")
for i in range(N1):
print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for i in range(N1):
print("{0:>7.2f} | ".format(Y1[i]), end="")
print("\n---------------------------------------------------------------------------------------------------------------------------------------------------------------------\nFFT3:")
#Y2 = np.absolute(Y2)
print("i = | ", end="")
for i in range(N1):
print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for i in range(N1):
print("{0:>7.2f} | ".format(Y3[i*4]), end="")
print("\n---------------------------------------------------------------------------------------------------------------------------------------------------------------------\nPer Hand:")
print("i = | ", end="")
Y = np.absolute(Y)
for i in k:
print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for i in k:
print("{0:>7.2f} | ".format(Y[i]), end="")
C:\Users\Fabia\AppData\Local\Temp\ipykernel_2544\3322552776.py:8: MatplotlibDeprecationWarning: The 'use_line_collection' parameter of stem() was deprecated in Matplotlib 3.6 and will be removed two minor releases later. If any parameter follows 'use_line_collection', they should be passed as keyword, not positionally. ax.stem(t3, y3, 'r',use_line_collection=True) C:\Users\Fabia\AppData\Local\Temp\ipykernel_2544\3322552776.py:19: MatplotlibDeprecationWarning: The 'use_line_collection' parameter of stem() was deprecated in Matplotlib 3.6 and will be removed two minor releases later. If any parameter follows 'use_line_collection', they should be passed as keyword, not positionally. ax.stem(f3, Y3, 'r',use_line_collection=True)
FFT1: i = | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | Y = | 250.00 | 226.53 | 163.32 | 79.55 | 0.00 | 53.15 | 67.65 | 45.06 | 0.00 | 45.06 | 67.65 | 53.15 | 0.00 | 79.55 | 163.32 | 226.53 | --------------------------------------------------------------------------------------------------------------------------------------------------------------------- FFT3: i = | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | Y = | 250.00 | 226.53 | 163.32 | 79.55 | 0.00 | 53.15 | 67.65 | 45.06 | 0.00 | 45.06 | 67.65 | 53.15 | 0.00 | 79.55 | 163.32 | 226.53 | --------------------------------------------------------------------------------------------------------------------------------------------------------------------- Per Hand: i = | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | Y = | 250.00 | 225.08 | 159.15 | 75.03 | 0.00 | 45.02 | 53.05 | 32.15 | 0.00 | 25.01 | 31.83 | 20.46 | 0.00 | 17.31 | 22.74 | 15.01 |
Das Signal wird wieder mit $f_{S4} = 16\text{ kHz}$ abgetastet, allerdings ist das Messfenster jetzt nur noch $60 \, T_S$ lang.
Signal und Zeitvektor erzeugen Sie am einfachsten durch Verkürzen der Signale aus dem letzten Abschnitt, also mit y4 = y3[:60] oder alternativ y4 = y3[:-4] (negative Indizes zählen vom Ende aus). Der Index, der das Ende markiert ("60" bzw. "-4") definiert dabei den ersten Wert, der nicht mit ausgegeben wird. Das ist z.B. bei arange(0,10,1) genauso und ist eine Folge des 0-based indexing (der erste Wert eines Arrays wird mit 0 angesprochen).
fig, ax = plt.subplots(**size)
N4 = 60
f_S4 = 16 * 1e3
t4 = np.arange(N4)/f_S4
y4 = y[:N4*40:40]
#FFT
Y4 = np.fft.fft(y4, N4)
Y4 = (np.absolute(Y4)/N4)
f4 = np.fft.fftfreq(N4, 1/f_S4)
#Plot FFT
ax.set_title("Y4")
ax.stem(f4, Y4, 'r',use_line_collection=True)
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{mV} \;\rightarrow$");
fig, ax = plt.subplots(**size)
ax.set_title("Y3")
ax.stem(f3, Y3, 'r',use_line_collection=True)
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{mV} \;\rightarrow$");
# ----------- Tabellenausgabe --------------
print("FFT1:")
print("i = | ", end="")
for i in range(N1):
print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for i in range(N1):
print("{0:>7.2f} | ".format(Y1[i]), end="")
print("\n---------------------------------------------------------------------------------------------------------------------------------------------------------------------\nFFT4:")
#Y4 = np.absolute(Y4)
print("i = | ", end="")
for i in range(N1):
print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for i in range(N1):
print("{0:>7.2f} | ".format(Y4[i]), end="")
print("\n---------------------------------------------------------------------------------------------------------------------------------------------------------------------\nPer Hand:")
print("i = | ", end="")
Y = np.absolute(Y)
for i in k:
print("{0:>7d} | ".format(i), end="")
print("\nY = | ", end="")
for i in k:
print("{0:>7.2f} | ".format(Y[i]), end="")
C:\Users\Fabia\AppData\Local\Temp\ipykernel_2544\933003921.py:15: MatplotlibDeprecationWarning: The 'use_line_collection' parameter of stem() was deprecated in Matplotlib 3.6 and will be removed two minor releases later. If any parameter follows 'use_line_collection', they should be passed as keyword, not positionally. ax.stem(f4, Y4, 'r',use_line_collection=True) C:\Users\Fabia\AppData\Local\Temp\ipykernel_2544\933003921.py:21: MatplotlibDeprecationWarning: The 'use_line_collection' parameter of stem() was deprecated in Matplotlib 3.6 and will be removed two minor releases later. If any parameter follows 'use_line_collection', they should be passed as keyword, not positionally. ax.stem(f3, Y3, 'r',use_line_collection=True)
FFT1: i = | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | Y = | 250.00 | 226.53 | 163.32 | 79.55 | 0.00 | 53.15 | 67.65 | 45.06 | 0.00 | 45.06 | 67.65 | 53.15 | 0.00 | 79.55 | 163.32 | 226.53 | --------------------------------------------------------------------------------------------------------------------------------------------------------------------- FFT4: i = | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | Y = | 266.67 | 18.52 | 26.52 | 62.62 | 212.93 | 55.77 | 51.29 | 113.09 | 99.64 | 34.91 | 28.87 | 81.28 | 16.67 | 4.41 | 1.45 | 0.00 | --------------------------------------------------------------------------------------------------------------------------------------------------------------------- Per Hand: i = | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | Y = | 250.00 | 225.08 | 159.15 | 75.03 | 0.00 | 45.02 | 53.05 | 32.15 | 0.00 | 25.01 | 31.83 | 20.46 | 0.00 | 17.31 | 22.74 | 15.01 |
In diesem Versuchsteil versuchen wir mit Hilfe der FFT Amplituden und Frequenzen von zwei Sinustönen zu schätzen mit $f_1 = 990\text{ Hz}$ und $A_1 = 100\text{ mV}$ sowie $f_2 = 1010\text{ Hz}$ und $A_2 = 2\text{ mV}$, überlagert ist eine gleichverteilte Rauschstörung im Bereich $\pm 1 \text{ mV}$.
Tipp: Diesen Versuchsteil können Sie auch mit pyfda durchführen, im Tab y[n] wählen Sie dazu ein Sinussignal mit überlagertem Rauschen aus und schauen sich das Spektrum im Untertab "Frequency" an. Enablen Sie "Stimulus X" und deaktivieren Sie "Response Y".
Die Anzahl der Datenpunkte $N$ für die FFT und die Abtastfrequenz $f_S$ können zunächst frei gewählt werden.
Welche Abtastfrequenz benötigen Sie mindestens?
Wählen Sie die kleinstmöglichen Werte für $N$ und $f_S$ so, dass kein Frequenzfehler und kein Leckeffekt auftritt. Tipp: Bestimmen Sie mit $T_{min} = L_1 T_1 = L_2 T_2$ zunächst das kürzeste Zeitfenster, in dem sowohl $f_1$ als auch $f_2$ periodisch sind. Danach können Sie mit $T_{Mess} = \frac{N}{f_S} = M T_{min}$ eine Bedingung für $N$ und $T_S$ ableiten
Wie groß ist das Verhältnis der Amplituden und wie groß das Verhältnis der Leistungen beider Sinustöne in dB?
Wählen Sie $f_S = 2500 \text{ Hz}$ und $N = 250$ für ein leichter ablesbares Spektrum.
Stellen Sie das Spektrum als (Amplituden-)Betragsspektrum im logarithmischen Maßstab dar, die Werte sollen relativ zum Maximum dargestellt werden. Zur Verbesserung der Darstellung können Sie mit np.fft.fftshift positive und negative Achse von Frequenzvektor und Spektrum vertauschen.
fig, ax = plt.subplots(**size)
N5 = 250
f_S5 = 2500
t = np.arange(N5)/f_S5
noi = 2 * np.random.rand(N5) - 1
y5 = 100 * np.sin(2 * np.pi * 990 * t) + 2 * np.sin(2 * np.pi * 1010 * t) + noi
#y5 = y5[:N5:40]
#FFT
Y5 = np.fft.fft(y5, N5)
Y5 = np.absolute(Y5)/N5
Y5 = 20*np.log10(Y5) #Umrechnen in dB
f5 = np.fft.fftfreq(N5, 1/f_S5)
#Y5 = np.fft.fftshift(Y5)
#f5 = np.fft.fftshift(f5)
#Plot FFT
ax.stem(f5, Y5, 'r',use_line_collection=True)
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
#ax.set_yscale('log')
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{dBmV} \;\rightarrow$");
C:\Users\Fabia\AppData\Local\Temp\ipykernel_2544\1150239569.py:20: MatplotlibDeprecationWarning: The 'use_line_collection' parameter of stem() was deprecated in Matplotlib 3.6 and will be removed two minor releases later. If any parameter follows 'use_line_collection', they should be passed as keyword, not positionally. ax.stem(f5, Y5, 'r',use_line_collection=True)
Die Abtastfrequenz ist jetzt fest eingestellt auf $f_S = 4 \text{ kHz}$, die Anzahl der Datenpunkte $N$ für die FFT soll eine Zweierpotenz sein, $N = 2^L$, für eine möglichst effiziente Berechnung der FFT.
Auch diesen Unterpunkt können Sie mit pyfda bearbeiten.
Der Leckeffekt soll jetzt durch Verwendung "weicherer" Fenster reduziert werden. Fensterfunktionen finden Sie unter scipy.signal.windows (https://docs.scipy.org/doc/scipy/reference/signal.windows.html). Vor Berechnung der FFT müssen Sie dazu Ihr Signal mit der Fensterfunktion multiplizieren, also z.B. y_win = y * sig.windows.hann(256, sym=False). Dabei wurde angenommen, dass auch Ihr Zeitsignal y eine Länge von $N = 256$ hat. Für die Spektralanalyse periodischer Signale (im Gegensatz zum Filterdesign und zur Analyse von impulsförmigen Signalen) muss sym = False gewählt werden.
Testen Sie den Einfluss der Fensterlänge $N = 256$, $N=512$ und $N=1024$ ... und den Effekt eines rechteckförmigen (boxcar bzw. gar keine Fensterung) und eines Hann-Fensters. Wie genau werden die Amplituden und Frequenzen geschätzt? Schauen Sie sich auch die anderen Fensterfunktionen an.
#Variablen
N6 = 1024
f_S6 = 4000
t = np.arange(N6)/f_S6
noi = 2 * np.random.rand(N6) - 1
y6 = 100 * np.sin(2 * np.pi * 990 * t) + 2 * np.sin(2 * np.pi * 1010 * t) + noi
#Windowing
y6_win = y6 * sig.windows.hann(N6, sym=False)
#y6_win = y6 * sig.windows.boxcar(N6, sym=False)
#y6_win = y6 * sig.windows.tukey(N6, sym=False)
#FFT
Y6 = np.fft.fft(y6, N6)
Y6_win = np.fft.fft(y6_win, N6)
#Nachverarbeitung
Y6 = np.absolute(Y6)/N6
Y6_win = np.absolute(Y6_win)/N6
#Y6 = 20*np.log10(Y6)
#Y6_win = 20*np.log10(Y6_win)
f6 = np.fft.fftfreq(N6, 1/f_S6)
#Y6 = np.fft.fftshift(Y6)
#f6 = np.fft.fftshift(f6)
#Plot FFT
fig, ax = plt.subplots(**size)
ax.set_title("Ohne Fenster")
ax.stem(f6, Y6, 'r', 'b', use_line_collection=True)
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{dBmV} \;\rightarrow$")
fig, ax = plt.subplots(**size)
ax.set_title("Mit Fenster")
ax.stem(f6, Y6_win, 'r', 'b', use_line_collection=True)
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y_1(F)\; / \;\mathrm{dBmV} \;\rightarrow$")
C:\Users\Fabia\AppData\Local\Temp\ipykernel_2544\953271478.py:35: MatplotlibDeprecationWarning: The 'use_line_collection' parameter of stem() was deprecated in Matplotlib 3.6 and will be removed two minor releases later. If any parameter follows 'use_line_collection', they should be passed as keyword, not positionally. ax.stem(f6, Y6, 'r', 'b', use_line_collection=True) C:\Users\Fabia\AppData\Local\Temp\ipykernel_2544\953271478.py:35: MatplotlibDeprecationWarning: Passing the markerfmt parameter positionally is deprecated since Matplotlib 3.5; the parameter will become keyword-only two minor releases later. ax.stem(f6, Y6, 'r', 'b', use_line_collection=True) C:\Users\Fabia\AppData\Local\Temp\ipykernel_2544\953271478.py:41: MatplotlibDeprecationWarning: The 'use_line_collection' parameter of stem() was deprecated in Matplotlib 3.6 and will be removed two minor releases later. If any parameter follows 'use_line_collection', they should be passed as keyword, not positionally. ax.stem(f6, Y6_win, 'r', 'b', use_line_collection=True) C:\Users\Fabia\AppData\Local\Temp\ipykernel_2544\953271478.py:41: MatplotlibDeprecationWarning: Passing the markerfmt parameter positionally is deprecated since Matplotlib 3.5; the parameter will become keyword-only two minor releases later. ax.stem(f6, Y6_win, 'r', 'b', use_line_collection=True)
Text(0, 0.5, '$Y_1(F)\\; / \\;\\mathrm{dBmV} \\;\\rightarrow$')
In diesem Versuchsteil schätzen wir das Fourierintegral des rechteckförmigen Impulses $y(t)$ mit $\Delta T = 250 \,\mu$s und $A = 1000\,\text{mV}$ (s. nächster Plot) mit Hilfe der FFT.
Bei welcher Frequenz hat die Fouriertransformierte die erste Nullstelle?
fig, ax = plt.subplots(**size)
Delta_T = 250e-6
A = 1000
t = np.arange(-200e-6, 200e-6, 1e-6)
ax.plot(t*1e6, A*np.where((t >= -Delta_T/2) & (t < Delta_T/2), 1, 0), 'r', lw=2)
ax.axvline(0, ls='-',color='k')
ax.set_xlabel(r"$t \; / \;\mu\mathrm{s} \;\rightarrow$")
ax.set_ylabel(r"$y(t) \; \mathrm{in \;mV} \;\rightarrow$");
Das zeitkontinuierliche Signal $y(t)$ hat eine endliche zeitliche Ausdehnung und endliche Energie, daher kann sein Amplitudenspektrum mit dem Fourierintegral beschrieben werden:
$$Y\left( f \right) = \int_{-\infty}^{\infty}y(t)\text{e}^{-\text{j} 2 \pi f t} \text {d}t$$Der Impuls ist symmetrisch zur y-Achse was die Berechnung deutlich vereinfacht:
$$ \begin{align} Y\left( f \right) &= A \int_{-\Delta T/2}^{\Delta T/2}\text{e}^{-\text{j} 2 \pi f t} \text {d}t = \left. \frac A{-\text{j} 2 \pi f} \text{e}^{-\text{j} 2 \pi f t}\right|_{-\Delta T/2} ^{\Delta T/2} = A\frac{\text{e}^{\text{j} \pi f \Delta T}- \text{e}^{-\text{j} \pi f \Delta T}}{2 \pi \text{j}f} = A\frac{\sin \pi f \Delta T}{\pi f} \\ &= A\Delta T \frac{\sin \pi f \Delta T}{\pi f \Delta T} = A\Delta T \text{sinc} f \Delta T \text{ mit } \text{sinc} f = \frac{\sin \pi f}{\pi f} \end{align}$$Das Ergebnis ist die Amplitudenspektrumsdichte in V/Hz.
Berechnet man die FFT anstatt des Fourierintegrals, überstreicht jeder Frequenzpunkt die Bandbreite $\Delta f = f_S/N$. Möchte man die FFT eines Energiesignals so skalieren dass man (in etwa) identische Zahlenwerte erhält wie bei der spektralen Amplitudendichte des ursprünglichen zeitkontinuierlichen Signals, so muss man mit $1/\Delta f = N/f_S$ multiplizieren. Da die FFT bereits mit $1/N$ skaliert ist, hebt sich der Faktor $N$ auf. Man muss also nur durch $f_S$ dividieren.
Für unendlich ausgedehnte Leistungssignale (periodische Signale, stationäre und nicht-stationäre Prozesse) konvergiert das Integral nicht, man nimmt stattdessen die Fourierreihe (periodische Signale) oder die spektrale Leistungsdichte, die über die Autokorrelationsfunktion berechnet wird (aber nicht hier und jetzt ...).
fig, ax = plt.subplots(**size)
Delta_T = 250e-6
A = 1000
f = np.arange(-2.2e4, 2.2e4, 1e2)
ax.plot(f, A * Delta_T*np.sinc(f * Delta_T), 'r', lw=2)
ax.axvline(0, ls='-',color='k'); ax.axhline(0, ls='-',color='k')
ax.set_title(r"Fourierintegral des Rechteckimpulses mit $\Delta T = 250\, \mu \mathrm{s}$")
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y(f) \; \mathrm{in \;mV/Hz} \;\rightarrow$");
Im folgenden soll die Fouriertransformierte des Impulses mit Hilfe einer FFT abgeschätzt werden, dazu wird der Impuls mit $L = 16$ Samples abgetastet.
Welche Abtastfrequenz $f_S$ ist dazu notwendig?
$f_S = $ 16 * 1/400 us = 40 kHz | T_mess = 400 us
Wie groß ist die Frequenzauflösung $\Delta f$? Wieviele Frequenzpunkte erhält man zwischen $f=0$ und der ersten Nullstelle $f_0$?
$\Delta f = $ f_s / N = 40kHz / 16 = 2.5 kHz -> f_0 = 4kHz -> f_0 / f_delta = 1.6 -> 1 (+1 wenn mit f=0)
Kann man die graphische Darstellung des Amplitudenspektrums zwischen $f= 0$ und der ersten Nullstelle $f_0$ durch Änderung der Abtastfrequenz verbessern?
Wie kann die graphische Darstellung des Betragsgangs verbessert werden ohne die Abtastfrequenz zu verändern?
Wie müssen die Ergebnisse skaliert werden, um physikalisch korrekte Werte für die spektrale Amplitudendichte in V/Hz zu bekommen (so wie theoretisch berechnet)? Tipp:
Berechnen Sie die FFT des Rechteckpulses mit $L=16$ und mit $N_{FFT} = 2^9 = 512$ Punkten.
Erstellen Sie einen Plot mit der FFT mit $L=16$ Punkten, $N_{FFT} = 512$ Punkten und dem idealen Amplitudenspektrum des Rechteckimpulses.
fig, ax = plt.subplots(**size)
f_S = 64e3
Delta_T = 250e-6
A = 1000
f = np.arange(-32e3, 32e3, 1e2)
f_16 = np.fft.fftshift(np.fft.fftfreq(16, d= 1/f_S))
f_512 = np.fft.fftshift(np.fft.fftfreq(512, d= 1/f_S))
y_id = A * Delta_T*np.sinc(f * Delta_T)
y_16 = A*np.ones(16)
Y_16 = np.fft.fftshift(np.abs(np.fft.fft(y_16))/f_S)
Y_512 = np.fft.fftshift(np.abs(np.fft.fft(y_16, 512))/f_S)
ax.plot(f, y_id, 'r', lw=2, label="$Y_{id}$")
ax.plot(f_16, Y_16, 'bo', label="$Y_{16}$")
ax.plot(f_512, Y_512, label="$Y_{512}$")
ax.legend()
ax.axvline(0, ls='-',color='k'); ax.axhline(0, ls='-',color='k')
ax.set_xlabel(r"$f \; / \;\mathrm{Hz} \;\rightarrow$")
ax.set_ylabel(r"$Y(f) \; \mathrm{in \;V/Hz} \;\rightarrow$");
(c) 2016 - 2021 Prof. Dr. Christian Münker
This jupyter notebook is part of a collection of notebooks on various topics of Digital Signal Processing. The latest version can be found at https://github.com/chipmuenk/dsp.
This notebook is provided as Open Educational Resource. Feel free to use it for your own purposes. The text is licensed under Creative Commons Attribution 4.0, the code of the IPython examples under the MIT license. Please attribute the work as follows: Christian Münker, Digital Signal Processing - Vorlesungsunterlagen mit Simulationsbeispielen, 2020.